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Abstract 

We adapt the simulated annealing algorithm to the search of periodic orbits for classical multi-electron 
atomic systems. This is done by minimizing the n-th return distance to the initial position on a Poincare 
surface of section under an energy constraint. Here we give evidence of the feasibility of the method by 
applying it to the helium atom in the ground state for one to three spatial dimensions. We examine the 
structure of the dynamics and connect its organization to the periodic orbits we have found. 

Keywords: Periodic orbits, Simulated annealing, Hamiltonian systems. 



1. Introduction 

Periodic orbits play a central role in the description and analysis of dynamical systems as they represent 
the skeleton of the dynamics It means that some important dynamical properties can be deduced 
from these orbits. By continuity of the flow a typical trajectory approaching a periodic orbit will mimic 
its dynamics. The time delay during which the trajectory is caught by the periodic orbit depends on 
the stability of the orbit: the more stable the periodic orbit, the longer the trajectory will stay in its 
neighborhood. As a result, the knowledge of the periodic structures and their stability properties enables 
one to predict the dynamical organization of the flow in their neighborhood. A fair amount of information 
is provided by the linear stability properties of these invariant structures. As an example, cycle expansions 
according to the length, stability or action of these orbits, are carried out to describe long time behavior 
such as averages of observables P, Q. For chaotic systems, the symbolic dynamics describes a hierarchy 
between periodic orbits, and the resulting expansion converges exponentially or superexponentially with the 
cycle length. For autonomous Hamiltonian systems, the eigenvalues of the Jacobian matrix from which the 
linear stability properties are determined come in quadruplet (A, 1/A, A*, 1/A*). In addition, there are at 
least two marginal eigenvalues corresponding to the time translation invariance (along the periodic orbit) 
and the energy conservation. These eigenvalues allow the classification of periodic orbits according to their 
linear stability properties. For instance, for Hamiltonian systems with two degrees of freedom, the periodic 
orbits can be sorted in three categories depending on their linear stability property: they are either elliptic 
(in general, linearly stable), hyperbolic (linearly unstable) or parabolic (linearly neutral). Elliptic periodic 
orbits are generally surrounded by an elliptic island inside which trajectories are trapped (and stay on 
invariant tori). In this region, the dynamics is mainly ruled by the central periodic orbit and the size 
of the island is determined by nonlinear stability properties. For hyperbolic periodic orbits, neighboring 
trajectories (linearly or locally) exponentially diverge in time. In general, these orbits are surrounded by 
a chaotic dynamics of stretching and folding since their stable and unstable manifolds intersect an infinite 
number of times to create a chaotic tangle. 

Of course the influence of various kinds of periodic orbits strongly depends on the problem at hand. For 
atomic and molecular systems, the typical duration of a process is short compared to the other physical 
processes. This is particularly true for systems driven by short and intense laser pulses. Periodic orbits 
longer than this typical duration (or even of the same order) will not influence drastically the dynamics. On 
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the contrary, periodic orbits much shorter than this duration will have a chance to trap the trajectory in its 
neighborhood and hence significantly affect the dynamical properties. 

Various algorithms have been developed for finding periodic orbits. Among them, some are commonly 
used in atomic systems with a few number of electrons : Newton- Raphson algorithm or modified version of it 
such as the damped Newton- Raphson algorithm offers high speed convergence provided one has a sufficiently 
accurate initial guess for the periodic orbit. The precision of this initial guess is a thorny problem since 
the basin of attraction shrinks exponentially with the instability and the length of the orbit 0, 0|- A 
solution to overcome this difficulty is to consider a parallelized version of the algorithm through a multi- 
shooting strategy [![. Independently of the chosen version, a crucial point in Newton- Raphson methods 
is the evaluation of the derivative of the trajectory which can be a delicate point from a numerical point 
of view. For a flow with dimension N, the tangent flow is given by the evolution of a N x N matrix 
which increases significantly the dimensionality of the problem at hand (even if in some cases, the structure 
of the problem may enable one to reduce the actual number of components). Some methods were set 
up to overcome this issue, and are generally based on relaxation algorithms 0, Overall, the goal is 

to set up a deterministic dissipative dynamical system which has the periodic orbit as its attractor (with 
hopefully a wide basin of attraction). Using the same philosophy, variational methods have been designed to 
determine periodic orbits. For instance, a variational principle combined with a Newton descent Q consists 
in setting up a (dissipative) fictitious-time dynamics in a space of loops such that it drives a loop to a true 
periodic orbit of the considered dynamical system. Similarly some other algorithms are designed based on 
minimizing a function (or a functional). For instance one can consider a Newton-Gauss method by looking 
at minimization of the distance between the starting and final points of an orbit on a particular surface 
intersecting the periodic orbit: Global vanishing minima correspond to periodic orbits. 

Independently of the chosen algorithm, the chance to see the algorithm converging strongly depends 
on the basin of attraction of the periodic orbit of the deterministic dynamical system. Often, there is a 
trade-off between the size of the basin of attraction with the speed of convergence of the algorithm. This 
is in particular the case for the Newton-Raphson and its damped version. In this article we propose a 
method which extends this basin of attraction by use of systematized trial and error converging procedure. 
In order to do that, we combine these deterministic algorithms with a Simulated Annealing (SA) algorithm 
to approximate the periodic orbits. In other words, the goal of this stochastic method is to enable the 
determination of initial guesses with sufficient accuracy to lay them into the basin of attraction of a fast 
converging algorithm like the Newton-Raphson algorithm. As a consequence of the underlying stochastic 
nature of the algorithm, it enables one to determine several different periodic orbits for the considered 
dynamical system by launching the algorithm several times. One of the advantages of the SA algorithm is 
that it does require neither the computation of the tangent flow nor a high accuracy in the integration of 
the trajectory. We show below that it makes this method a tool of choice for systematic search for periodic 
orbits in phase space. 

To illustrate the feasibility of the SA algorithm for finding periodic orbits, we consider the classical 
helium atom in its ground state. The dynamics is modeled by the following Hamiltonian with soft Coulomb 
potentials (loj : 

ut \ IPil 2 Jp 2 | 2 2 2 I 

W(xi,x 2) Pi,P2) = — — + — ; 1 + i =■ (!) 

Xl | 2 + a 2 V|x 2 | 2 + a 2 V /|x 1 -x 2 | 2 + 62 

where Xj is the position of the z-th electron (the nucleus being set at the origin), and Pi is its canonically 
conjugate momentum. The norm • is the Euclidian norm. Depending on the dimensionality of the model, 
the vectors Xj and Pi can be considered in K, R 2 or R 3 . Tuning the softening parameters a and b enables 
to reproduce various atomic configurations. The parameter a is chosen as to reproduce the ground state 
energy (defined as the sum of the first and second ionization potentials) while preventing any self ionization 
of the atom. The parameter b is chosen as to allow significant energy exchange between the electrons when 
they are close to each other. For helium, one usually considers a = b = 1 and the ground state energy is 
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The dynamical organization corresponding to Hamiltonian ([T]) is crucial for understanding the (multiple) 
ionization processes of these atoms driven by an intense and short laser pulse. In Ref. [12j | it was argued 
from a one dimensional model that the dynamics of Hamiltonian (fTJ is mainly organized by a reduced set of 
periodic orbits which naturally places one electron close to the nucleus and one further away. This distinction 
between the two electrons is crucial when understanding the action of an intense linearly polarized laser 
field when it is applied to the system : The electron far away from the nucleus is picked up by the field 
and hence quickly ionizes while the other electron remains bound to the nucleus. This distinction between 
the two electrons is at the origin of the classical interpretation of the recollision or the " 'three-step" ' model 
which is the keystone of strong- field physics [l3T - [l5j ]. The question we are addressing here is whether this 
emerging picture of an inner (close to the nucleus) and an outer electron (further away) still exists in two 
or three dimensions, or if it is just an artifact of the one dimensional model which is very peculiar given the 
presence of the nucleus. We show below that this dynamical picture still persists and and it does so because 
of the organization of short periodic orbits in phase space. 

In Sec. [21 we describe the algorithm we use to determine periodic orbits. In Sec. 03 we analyze the lay-out 
of these periodic orbits in phase space and their implication for the dynamical organization of the helium 
atom in the ground state. 

2. Simulated annealing algorithm to determine periodic orbits 

The Simulated Annealing (SA) algorithm [l6| is a metaheuristic algorithm designed for optimization 
under constraints which generalizes the ideas developed from the Metropolis algorithm [l]} ■ The main idea 
of the algorithm is to automate a trial and error search for a better solution within a controlled neighborhood 
whose diameter is successively reduced. 

Given a dynamical system X = F (X), we first consider a Poincare section E and the associated Poincare 
map $ : E M> E. If a trajectory cuts E in a finite set of points (Xi, X2, . . . , X„) such that $ (X„) = Xi 
and $(XJ = X, + i, 1 < i < n it corresponds to a periodic orbit with period n on the section (in other 
words $™(Xi) = Xi). The periodic orbit search is obtained by minimizing the distance between the 
starting point and the n-th return on E using the SA algorithm. From a random initial guess X(°) G E we 
build up a sequence {X( fc '} fc6N G E N which aims at minimizing the Euclidean distance between X( fc ) and 
$ n (X( fe )). If the SA algorithm converges, the limit X^ 00 ) of the sequence exists and fulfills the condition 
$ n (X(°°)) = X^ 00 ', therefore corresponding to a periodic orbit (with period n on E). For the problem at 
hand, finding periodic orbits of Hamiltonian (TTJ) in the ground state, an additional constraint to take care 
of during the SA algorithm is the energy of the system (which is constant for autonomous Hamiltonian 
systems). There are several ways of doing it : The first one is to ensure that the points X^ remain on the 
energy surface after each step of the SA algorithm by an additional procedure which adjusts one coordinate. 
The second way is a specific prescription on the function to minimize which only ensures that the limit 

x(°°). 

For a given function / : T — > R, where T is a subset of E, called cost function, the goal is to approximate 
the global minimum of / over T. First a random initial condition X^ ' G T is chosen to initiate the SA 
algorithm and we define two parameters To, the initial temperature, and k G M + which serves as error 
tolerance. 

The SA algorithm proceeds in successive stages of annealing and cooling steps. 

• During the annealing process, a random initial guess X^ is chosen in the neighborhood of X^ ) 
in r (the size of the neighborhood is controlled by To; for simplicity one usually chooses a ball with 
diameter To). It corresponds to a melting of the system with a controlled temperature Tq. If / (X*- 1 ') < 
/ (X^ -*) then X^ 1 ) replaces X^ ^ in the following step. If not, X^ replaces X^ ^ with a probability 
exp(-(/(X«) - /(X(°)))/kT ). Otherwise, the same X^ ) is kept in the next step. The process of 
melting for a fixed temperature is iterated until the system reaches a steady state. 

• The cooling consists in decreasing the temperature to a new temperature T± < To, realized at the 
end of the annealing process. Then a new iteration of annealing with the updated temperature T\, is 
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performed. 



The successions of annealing and cooling steps are iterated until the system reaches a global steady state 
which is the resulting approximation of the global minimum. Since the size of the neighborhood where the 
perturbation is performed is controlled by the temperature, it has to shrink with iterations of cooling. At 
the limit for zero temperature, the neighborhood of control is reduced to a single point. 

The constant k is used to avoid the algorithm to be trapped in local minima. It has to be chosen carefully 
as two critical situations may arise: 

• If k is too large, the algorithm will not converge since too many "bad" perturbations are allowed. 

• If k is too small, one runs the risk to be trapped in local minima. 

The cooling law also has to be chosen carefully to allow the algorithm to sufficiently explore the region around 
the initial guess. If the cooling is too fast, the system is frozen in its current position and the algorithm may 
not converge. If the cooling is too slow, the algorithm spends a lot of time in useless explorations which 
slows down the process. Usually, an exponential decrease is considered for the cooling by multiplying the 
current temperature with a constant, i.e. T n+ \ = aT n where < a < 1. 

Recall that our goal is to find periodic orbits for Hamiltonian ([1]) with the constraint to be on the ground 
state, and to do that, we are looking for n periodic points on a Poincare surface. For instance, we choose the 
Poincare surface E of equation x\ — with i,\ > 0, where x% is the first component of x in the canonical basis 
(other surfaces which intersect the ground state energy surface might also be suitable). For convenience, 
we denote X = (xi, x 2 , pi, P2) the vector position in phase space. Our problem is to find a periodic point 
X(°°J on the Poincare (and the ground state energ y) surface such that X^ = (X(°°)). We define the 
cost function as the n-th return distance function 



/ : X |X - $ n (X) 



(2) 



such that periodic orbits with period n correspond to positions where the global minimum of the function 
/ is reached and thus may be investigated using the SA algorithm. This function is defined on T which is 
the intersection of the ground state energy surface {X s.t. H (X) = E g } with the Poincare section E. 
Schematically, the algorithm is written as 



°/ Initialization 
Xold 

DistError = 
Coollter = 



InitCondO 

ReturnDist(Xold) 





7. Random initial condition generation 

7o n return distance from initial condition 

7» Cooling iteration counter 



% Simulated annealing algorithm 

while (Coollter < #CoolIter) and (DistError < #DistError) 

Meltlter = °L Melting iteration counter 

while (Meltlter < #MeltIter) and (DistError < #DistError) 

Xnew = Xold + Tprt*2* (randO - . 5) °L (*) random perturbation of current guess 
Dist = ReturnDist (Xnew) 7o n return distance for perturbed condition 



if Dist < DistError 

Xold = Xnew 

DistError = Dist 



7o Improved guess 



7o Allowed not improving perturbation 

else if (exp(-(Dist-DistErr)/(#ErrTol*Tprt)) < randO) 
Xold = Xnew 

DistError = Dist 
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end if 



Meltlter = Meltlter + 1 
end while 

Tprt = Tprt*#TprtFact 
Coollter = Coollter + 1 

end while 

where user-defined constants are labeled by # and % starts comments. The function ReturnDist plays the 
role of the function / as defined by Eq. ©. Finally, the function rand is a random generator uniformly 
distributed over [0, 1]. The melting process (*) is draffcd in its simplest way above and it has to be followed 
by a projection on T, which has to be chosen accordingly to the problem at hand. In our case, after the 
melting (*), the new initial condition is neither on T, nor on the ground state energy. The projection on E 
is simply done by setting x% = 0. The projection on the ground state is done by changing p Xi i accordingly 
to the energy constraint and we select the positive solution for p Xt i in order to be sure to remain on the 
Poincare section. We also define a threshold distance under which the algorithm is considered as converged. 

The projection of the perturbed position on the ground state energy surface ensures that the series of 
guesses stay on the energy surface at each iteration. However, one can also see the ground state energy as 
an additional constraint to minimize together with the n-return distance. In this case, one adds the energy 
error to the cost function 

/:X^|X-$"(X)|+ 7 |H(X)-£ 3 |, (3) 

where 7 > is a constant. Then, after each perturbation, we only project the new guess on the Poincare 
surface S. We have successfully implemented both methods for n £ {1, 2, 3}. All numerical results displayed 
in this paper and the values given for the parameters correspond to the projection (on both the Poincare 
and ground state energy surfaces) method. 

Many improvements of the algorithm can be implemented depending on the problem at hand. For 
instance, there is no requirement for choosing the initial condition randomly. If one has an approximate 
guess for the location of a periodic orbit, or any subset of phase space where it is included (e.g., given the 
symmetries of the problem), it is more suitable to start the SA algorithm with this particular point. In such 
a case, it can be interesting to set the initial temperature cooler than for a random initial guess. 

For improving the convergence, it is more efficient to combine the SA algorithm with a second method 
based on a Newton-Raphson method for instance. Here we use the function f solve of Matlab, which is a 
Newton-Gauss minimizing scheme. 

3. Application to two electron atoms with soft Coulomb potentials 

For two spatial dimensions the complexity of the system is increased compared to one dimensional models 
since the number of degrees of freedom is increased from 2 to 4. In addition, for two spatial dimensions, 
the system has a continuous symmetry which corresponds to a global invariance by rotation of the atom in 
phase space. The associated conserved quantity is the total angular momentum xi x pi +X2 x P2, where x 
is the cross product in R 2 : x x p = xp y — yp x - It means that any periodic orbit generates an infinite family 
of periodic orbits deduced from this single one by a rotation of the atom in phase space. In addition, the 
system also has some discrete symmetries (the exchange of the role of the two electrons for instance). To 
simplify the analysis, in what follows, we identify as a single one all periodic orbits which can be deduced 
by one of those symmetries from another periodic orbit. We notice that these periodic orbits have the 
same period and the same linear stability properties. In the numerical implementation, we consider that a 
periodic orbit has converged when the distance in phase space between the starting and final points on the 
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Poincare section is smaller than 10~ 10 . We consider a periodic orbit to be a new one when the difference in 
period or in the spectrum of the monodromy matrix is larger than 10 -2 from any already known periodic 
orbit in the list. 

For actual use of the SA algorithm we set the initial temperature around 5. A tolerance parameter 
k = + (a perturbation is accepted if an only if it reduces the cost function) is enough for our problem 
when we perform the projection on the ground state energy surface. It should be noted that if we add the 
energy error to the cost function [see Eq. ©], k«5x 10~ 2 is a better choice. For the cooling we set the 
temperature ratio to a — 0.75. Finally, we iterate the melting step a fixed number of times (3000 melting 
steps during the annealing process for each temperature). We stop the cooling when the temperature has 
reached a threshold value (here the threshold temperature is 10 -5 ). We also define a threshold d cr it for 
which the algorithm is stopped as soon as the current guess X satisfies / (X) < d cr a, and here d cr it = 10 -3 . 

For actual use of the SA, there are two main strategies when choosing the temperature ratio and the 
number of melting steps. One strategy is to choose a fast cooling law (i.e. a small) and a large number of 
melting steps for each annealing. Another strategy is to consider a slow cooling rate (i.e. a ~ 0.95) and a 
much reduced number of melting per annealing iteration. For the problem at hand, we have used the first 
strategy when considering a projection on the ground state energy surface, and we have used the second 
strategy when adding the energy error to the cost function. 

The total angular momentum is a conserved quantity of the dynamics and to model a configuration closer 
to the quantum ground state, one may wish to impose it to be zero. That becomes a new constraint for 
our system and it decreases the actual dimensionality of the problem at hand, but still it is worth looking 
for periodic orbits which may organize this particular subclass of the dynamics. To do so, there are two 
straightforward ways to adapt the algorithm presented here. The first one consists in adding the total 
angular momentum to the cost function (in the same way as for the energy constraint) 

/ : X .->• |X- $" (X)| + 7 |xi x pi + x 2 x p 2 | , 

where 7 is constant. The total angular momentum is then considered as an additional constraint and the 
SA will try to minimize it together with the n-th return distance function. Another solution is to project, 
after each melting iteration, the guess on the intersection of the ground state energy and zero total angular 
momentum surfaces. We have successfully adapted our algorithm to find one dimensional, three dimensional 
and zero total angular momentum two dimensional (considering the projection alternative) periodic orbits. 

Since we are interested in finding as many periodic orbits as possible up to a finite period due to the 
finite duration of laser pulses in experiments, we launch a large number of times the algorithm, to explore 
the whole accessible phase space. We do it for periodic orbits of Hamiltonian ((TJ) which have up to three 
intersections on the Poincare section. Because of the rotational invariance of the system, the monodromy 
matrix associated with a periodic orbit has an additional pair of eigenvalues 1 (in addition to the pair 
associated with time-translation along the orbit). As a result, for two spatial dimensions, four eigenvalues 
are equal to one in the monodromy matrix and the four others determine the linear stability properties. 
Consequently the system can exhibit nine kinds of linear stability, (elliptic, parabolic or hyperbolic)- (elliptic, 
parabolic or hyperbolic). However, for numerical computations, is it difficult to distinguish precisely the 
parabolic features, e.g., between a parabolic and an elliptic or a weakly hyperbolic periodic orbit. When 
checking the linear stability of the periodic orbits we have found, we arbitrarily group together parabolic 
and elliptic features (it is worth noting that we have found no fully parabolic periodic orbit). 

For the case of two spatial dimensions for Hamiltonian ((TJ, we have found a total of 155 periodic orbits 
with period smaller than 65 a.u. whose repartition is the following one : 43 periodic orbits with period one, 
74 periodic orbits with period two (which cannot be reduced to period one, that is to say, prime periodic 
orbits) and 38 periodic orbits with period three. Regarding their stability properties, 22 are elliptic-elliptic 
or parabolic-elliptic, 71 are elliptic-hyperbolic or parabolic-hyperbolic and 62 are hyperbolic-hyperbolic. In 
Fig. [TJ we display the period T as a function of the modulus of the largest eigenvalue A max of the monodromy 
matrix for each periodic orbit determined by the algorithm. By no means are these periodic orbits the only 
ones in phase space. However, these are the ones obtained by launching a large number of times the SA 
algorithm (with the restrictions mentioned above). 
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Figure 1: Period T versus the norm of the largest eigenvalue A max of the monodromy matrix for the detected periodic for 
Hamiltonian {TJ with two spatial dimensions. Elliptic-elliptic and parabolic-elliptic periodic orbits are labeled with red crosses, 
elliptic-hyperbolic and parabolic-hyperbolic ones with cyan squares and hyperbolic-hyperbolic ones with blue circles. Hyperbolic 
periodic orbits into the dashed black rectangle are depicted in Fig. [2] 



There are five periodic orbits which are weakly hyperbolic and with relatively low period, and which 
are likely to organize the dynamics of Hamiltonian ([T]) (see dashed black rectangle in Fig. [T|). This is so 
because all of those periodic orbits have a small period (around 15 a.u.) which is consistent with the short 
duration of laser pulses (typically of a few hundreds of a.u. which correspond to few tens of fs) and are 
weakly hyperbolic. These five hyperbolic periodic orbits are well separated, in the (A max , T) diagram, from 
the other periodic orbits which are significantly longer and more hyperbolic. We display these five periodic 
orbits in the configuration space (x, y) in Fig. [5] We notice that all of them are composed of one electron 
close to the nucleus and another further away. A closer inspection at Fig. [1] also shows a group of four 
weakly hyperbolic periodic orbits with a larger period (around 30 a.u.). Looking at their shape in phase 
space (not shown here) they are also composed of one electron close to the nucleus and the other one further 
away. It means that for a typical trajectory of Hamiltonian (TT]) under the influence of these periodic orbits 
it is possible to identify, at any time, an inner (close to the nucleus) and an outer (further away) electron 
with possible exchanges of the roles of the two electrons as the trajectory visits different areas of influence. 
The observation that all small period and weakly hyperbolic periodic orbits are clearly composed of an inner 
and an outer electrons, supports what we have shown for the model with one spatial dimension where we 
identified four dominant periodic orbits in the organization of phase space (which actually reduce to only 
one orbit by symmetry) [l2l Il8j|. However, for models with two spatial dimensions, the variety of periodic 
orbits and the invariance by rotation make it difficult to identify a reduced number of dominant orbits in 
the skeleton of the dynamics. In order to show that a typical trajectory has an inner and an outer electron, 
we consider two distances in phase space as functions of time : one from the nucleus and the other one 
from the boundary of the admissible phase shape for a typical trajectory of Hamiltonian (TTJ) as shown in 
Fig. [3] At any time an inner and an outer electron are identified with quick and frequent exchanges of their 
roles : The inner electron is close to the nucleus and thus have a short distance from it while the outer one, 
further away will be closer to the edge of the accessible phase space. It is worth noticing that due to the 
exchanges of the roles of the two electrons, the invariance by rotation and the dimensionality of the problem, 
a projection on the configuration plane does not clearly show any significant organization of the dynamics, 
contrary to the case with one spatial dimension which clearly shows the specific role played by mainly one 
periodic orbit (and its symmetric orbits). 

The same study can be carried out in three spatial dimensions. The numerical results suggest that in 
most of the periodic orbits, the two electrons are co-planar or nearly co-planar (but not all). One of them is 
depicted in Fig. [4l These periodic orbits display a similar organization with one electron close to the nucleus 
and one further away. As a consequence, from a systematic inspection of the periodic orbits in phase space, 
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Figure 2: Projection in the two dimensional configuration plane (x,y) of the five hyperbolic periodic orbits inside the dashed 
black rectangle in Fig. [T] Periodic orbits are sorted in increasing norm of the largest eigenvalue of the monodromy matrix, 
from left to right and up to down. Each electron is labeled with a different color (blue or red). 

we are able to draw a picture of the dynamics consisting of one electron close to the nucleus and another 
one further away (with fast exchanges between the two) , confirming a result obtained with a reduced system 
with only one spatial dimension. 
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